clc; clear; close all;
% 在四个不同的 k 值下，变量 x 和 y 的迭代序列

% 参数
mu = 3.5;
x0 = 0.1;
y0 = 0;
z0 = -1;
k_vals = [0.02, 0.15, 0.16, 0.2];  % 四个典型 k

N_transient = 500;  % 瞬态迭代步数
N_plot = 30;        % 用于绘制的迭代步数

figure('Position',[100 100 900 600]);

for i = 1:4
    k = k_vals(i);
    
    % 1) 瞬态迭代
    state = [x0, y0, z0];
    for n = 1:N_transient
        [dx, dy, dz] = mclm(state, mu, k);
        state = [dx, dy, dz];
    end
    
    % 2) 记录接下来 N_plot 步迭代的 (x, y)
    seqX = zeros(1, N_plot);
    seqY = zeros(1, N_plot);
    for n = 1:N_plot
        [dx, dy, dz] = mclm(state, mu, k);
        state = [dx, dy, dz];
        seqX(n) = dx;
        seqY(n) = dy;
    end
    
    % 3) 构造横坐标，使其从 501 ~ 530
    xAxis = (N_transient + 1) : (N_transient + N_plot);
    % xAxis(1) = 501, xAxis(end) = 530
    
    % 4) 绘制子图
    subplot(2,2,i);
    plot(xAxis, seqX, 'r-', 'LineWidth',1.2); hold on;
    plot(xAxis, seqY, 'b--','LineWidth',1.2);
    xlabel('n');
    ylabel('x_n, y_n');
    
    legend('x','y','Location','best');
    
    % 子图标题
    switch i
        case 1
            title('(a) k=0.02');
        case 2
            title('(b) k=0.15');
        case 3
            title('(c) k=0.16');
        case 4
            title('(d) k=0.20');
    end
    grid on;
end
